Contrasting assembly mechanisms and drivers of soil rare and abundant bacterial communities in 22-year continuous and non-continuous cropping systems

Despite many studies on the influence of cropping practices on soil microbial community structure, little is known about ecological patterns of rare and abundant microbial communities in response to different tobacco cropping systems. Here, using the high-throughput sequencing technique, we investigated the impacts of two different cropping systems on soil biochemical properties and the microbial community composition of abundant and rare taxa and its driving factors in continuous and rotational tobacco cropping systems in the mountain lands of Yunnan, China. Our results showed that distinct co-occurrence patterns and driving forces for abundant and rare taxa across the different cropping systems. The abundant taxa were mainly constrained by stochastic processes in both cropping systems. In contrast, rare taxa in continuous cropping fields were mainly influenced by environmental perturbation (cropping practice), while governed by deterministic processes under rotational cropping. The α-diversity indices of rare taxa tended to be higher than those of the abundant ones in the two cropping systems. Furthermore, the network topologies of rare taxa were more complex than those of the abundant taxa in the two cropping systems. These results highlight that rare taxa rather than abundant ones play important roles in maintaining ecosystem diversity and sustaining the stability of ecosystem functions, especially in continuous cropping systems.

Recent studies have shown that rare taxa may counteract environmental disturbances by increasing the functional redundancy of soil microbiota under environmental stress 13,14 , which remains to be empirically confirmed in different farming systems.
Rare taxa have a different response to environmental changes from abundant taxa 14,15 . Studies have observed that there were divergent assembly processes for abundant and rare taxa in oil-contaminated soils 15 and in agricultural soils 16 . Deterministic processes refer to biotic and abiotic factors determine microbial communities, while stochastic processes regard that microbial communities are mainly governed by random changes 17,18 . Many studies show that environmental factors alter the relative influences of stochastic and deterministic assembly processes [19][20][21][22] . For example, abundant taxa have been found to be constrained by stochastic processes in the surface layer of the northwestern Pacific Ocean 23 , while rare taxa have been found to be governed by stochastic in oil-contaminated soils 15 . At present, it is still unclear the primary forces that influence the assembly of rare and abundant taxa in continuous tobacco cropping soils.
The aims of this study were (1) to determine whether abundant and rare taxa show different responses to crop rotation and continuous cropping practices, (2) to investigate which taxa play a major role in the resilience to the continuous cropping, i.e., abundant or rare taxa, (3) to explore the relative importance of environmental variables influencing the assembly of rare and abundant sub-communities in the crop rotation and continuous cropping systems. For this, based on high-throughput Illumina paired-end sequencing technologies we analyzed the bacterial small-subunit ribosomal RNA (16S rRNA) gene sequences to determine both rare and abundant bacterial lineages in continuous and rotational tobacco cropping fields in the Yunnan Province, China. We hypothesized that: (1) the abundant and rare taxa exhibit distinct responses and assemblage patterns in continuous and noncontinuous tobacco soils; (2) the co-occurrence patterns of abundant and rare taxa would show divergence in both continuous and non-continuous tobacco soils, and their interactions will be weakened at the continuous tobacco fields due to continuous-cropping obstacle; (3) the underlying factors and driving forces regulating community assemblages vary across the continuous and non-continuous tobacco fields.
General distribution patterns of rare and abundant taxa. After quality filtering and the removal of chimeric sequences, the Illumina V4-V5 16S rRNA data set contained 1,460,302 quality sequences (range 124,400-477,309 sequence reads per sample) and 17,165 OTUs (Table S1). In both continuous cropping and crop rotation soils, the abundant taxa constituted a very low proportion of OTUs (mean 3.28% and 3.40%), but accounted for 41.30% and 41.81% of the average relative abundance in each sample. The rare taxa constituted a high proportion of the OTUs (mean 68.44% and 67.39% OTUs), while they contributed to an average of only 19.25% and 18.02% of the relative abundance in each sample (Table S1). The Goods coverage was 0.86 to 0.99, indicating that the number of sequence reads was sufficient to capture most taxa in each sample (Table 1). α-diversity (both the richness (Chao 1) and evenness (Shannon index) of the bacterial community in crop rotation soils tended to be higher than those in continuous cropping soils. The rare taxa showed higher α-diversity, both in terms of Chao1, Observed richness, and Shannon than the abundant taxa in the two cropping systems ( Table 2). The OUT richness of rare taxa was significantly higher than that of the abundant taxa in both continuous cropping and crop rotation soils (Fig. 1A,B). The evenness of rare taxa and abundant taxa was significantly higher than that of the total taxa in both continuous cropping and crop rotation soils, and the evenness of rare taxa was the highest (Fig. 1C,D). In the continuous cropping, the top 5 most dominant phyla were Actinobacteria, Acidobacteria_Gp4, Alphaproteobacteria, Acidobacteria_Gp6, and Unassigned. The top 5 most dominant phyla in the crop rotation were Actinobacteria, Acidobacteria_Gp4, Acidobacteria_Gp6, Sphingobacteriia, and   www.nature.com/scientificreports/ Unassigned. In both continuous cropping and crop rotation, the abundant taxa were more frequently Actinobacteria, Acidobacteria_Gp4, Acidobacteria_Gp6, Betaproteobacteria, Gemmatimonadetes, and Anaerolineae. The rare taxa were Alphaproteobacteria, Unassigned, Deltaproteobacteria, Planctomycetia, Gammaproteobacteria, and Chloroflexia. Acidobacteria_Gp3 and Chthonomonadetes were dominant in the rare sub-community in continuous cropping soils, but there was no significant difference in them between the rare and the abundant taxa in the rotation cropping soils. Bacilli dominated in the abundant taxa in continuous cropping soils, but it dominated in the rare taxa in crop rotation soils. There was no significant difference in Sphingobacteriia in the abundant and rare taxa in continuous cropping soils, but it dominated in the abundant and rare taxa in the crop rotation soils (Fig. 2). Different cropping modes had a significant influence on the abundance of total, abundant and rare taxa ( Fig. 3A-C). The results of NMDS showed that the similarity of the total, abundant and rare taxa between continuous cropping and crop rotation was low ( Fig. 3D-F). Rare taxa had higher community dissimilarity than abundant taxa in both cropping modes (Fig. 3F). Null model analysis showed that the influences of determinism (i.e. βNTI ≥ 2 or βNTI ≤ − 2) and stochasticity (i.e. − 2 < βNTI < 2) on assembly varied between the continuous cropping and crop rotation systems. Most of the β-NTI values of total and the rare taxa in the continuous cropping and crop rotation were higher than 2 (Fig. 4), while the β-NTI values of the abundant taxa were between− 2 and 2, indicating that the deterministic processes dominate the communities of the total and rare taxa in continuous cropping soils, while the random processes dominant in abundant taxa. In addition, the β-NTI values of the three sub-communities were between − 2 and 2 in crop rotation soils, indicating the three sub-communities were governed by the stochastic assembly in rotation soils. The β-NTI values of the total and rare taxa were higher than those of the abundant taxa (Fig. 4).
Co-occurrence patterns of the three sub-communities. The topological characteristics of the cooccurrence network showed that the nodes, modularity, clustering coefficient, average path length, network diameter, and average degree were all higher in the crop rotation soils than in the continuous cropping, while continuous cropping had more edges (Edges) ( Table 3). Compared with those in continuous cropping, the rare taxa in the crop rotation soils had more nodes and edges as well as higher modularity, average path length, and network diameter ( Table 3). The number of nodes, edges, modularity, and cluster coefficient of rare taxa was all higher, while the average path length, network diameter, and average degree were all lower than those of abundant taxa in both continuous cropping and crop rotation soils (Table 3). In addition, in Empirical networks, compared with the abundant taxa network, the rare taxa network was broader and had higher number of modules, modularity, and average geodesic distance (GD), while the abundant taxa had higher average connectivity (avgK) and average clustering coefficient (avgCC). In Random networks, GD and modularity in rare taxa networks were higher than those in the abundant taxa (Table 4).
In the two cropping modes, the nodes in the network graph mainly belong to 10 phyla: Acidobacteria, Gemmatimonadetes, Actinobacteria, Bacteroidetes, Chloroflexi, Plantctomycetes, and Proteobacteria (Fig. 5). The total taxa were more complex, while the abundant and rare taxa were more discrete in cropping rotation soils than those in continuous cropping soils. Compared with the abundant taxa, the rare taxa were more complex, as their structure was wider, and they had more nodes and edges and modularity in both continuous cropping and crop rotation soils (Fig. 5).
In the co-occurrence network of soil bacteria, most of the OTUs were divided into peripheral nodes (Peripherals) (Fig. S1). In continuous cropping soils, module hubs mainly belong to Acidobacteria_Gp4, Acidobacteria, Gemmatimonadales, Acidobacteria_Gp6, Deltaproteobacteria, Betaproteobacteria, Sphingobacteria, and Alphaproteobacteria, and most of these nodes had high intra-module connectivity (Zi) but the connectivity between modules (Pi) was low. In addition, connectors in the continuous cropping mainly belong to Actinobacteria, Acidobactena_Gp4, candidate_division, and Chloroflexi. These nodes had high inter-module connectivity (Pi) and low intra-module connectivity (Zi) (Fig. S1).
Relationship between microbial sub-communities and soil properties. The first two axes of the RDA analysis accounted for 29.66% and 11.39% of the variance in the total taxa, 48.31% and 13.01% in the abundant taxa, and 26.14% and 11.34% in the rare taxa. Among all the soil properties involved in the RDA analysis, AP, TOC, SOM, TP, TN, and TK were closely related to the soil bacterial communities in the two cropping patterns (Fig. 6). The Mantel test results indicated that the total microbial taxa were correlated with soil pH, TOC, SOM, TN, TP, TK, NH 4    www.nature.com/scientificreports/ while crop rotation could improve soil properties and soil microbial community 26,27 . We also found that the soil properties and microbial parameters were greater in crop rotation soil than the continuous cropping soil. As the soil TN, TP, TK, and all available nutrients (NH 4 + -N, NO 3 --N, AP, and AK) were obviously higher in crop rotation systems, which may relate to the change of soil pH under different cropping systems. Acidification could reduce the availability of nutrients and thereby reduce soil organism activity and biomass in soil [28][29][30] . Moreover, the improvement of soil texture and increased soil microorganisms are also tightly related to the enrichment in soil nutrient availability 31,32 . This indicates that soil nutrient contents were improved under crop rotation systems.
We found that the bacterial diversity was also greater in crop rotation soil than continuous cropping soil, and such responses may be linked to the soil physicochemical properties (e.g. soil pH, moisture, and nutrients content) 33,34 . As environmental variation is a key decisive factor of the distribution and abundance of microbes 35 . The relatively lower soil pH may decrease the bacterial diversity under continuous cropping systems, as soil acidity  www.nature.com/scientificreports/  www.nature.com/scientificreports/ could change the stability of bacterial cell membranes and thus inhibits bacterial growth 36,37 . In addition, the high content of soil organic C in crop rotation systems, to some extent, likely influences bacterial diversity, activities, and community compositions. A meta-analysis from Venter et al. 38 also revealed that increased microbial richness and diversity with crop rotation. In our study, soil microbial communities were further classified into rare and abundant microbial taxa. In general, abundant taxa had a higher capacity in competing for nutrient resources and thus not sensitive to environmental changes, and the opposite is true for rare taxa 39,40 . However, our results showed that the rare sub-communities have higher richness and α-diversity than the abundant taxa in both continuous cropping and crop rotation systems, indicating that rare taxa were, to some extent, not sensitive to the different cropping patterns, compared with abundant taxa 13,41 . The reason may be that the rare taxa may contribute to the resilience of the microbial community and sustain community function as an insurance source in continuous cropping soils due to their high functional redundancy and resistance to stress 14,40 . Especially, continuous cropping has strong negative effects on abundant taxa. Furthermore, greater community dissimilarity across the cropping systems was detected in the rare biosphere than in the abundant and total taxa, indicating the rare taxa may have a flexible taxonomic profile 11 . Actinobacteria, Acidobacteria_Gp4, Alphaproteobacteria, Acidobacteria_Gp6 were the predominant bacterial phyla in all the soil samples, which are known to be among the dominating phyla in different agricultural  www.nature.com/scientificreports/ practices 42 . Previous studies suggest that Actinobacteria is beneficial bacteria for plant growth 43 . Some research showed that the abundance of Actinobacteria decreased in continuous cropping soils 44 but increased in a crop rotation system 45 . Proteobacteria and Actinobacteria are copiotrophic in nature, which is known to have high nutritional requirements 46,47 , whereas Acidobacteria belongs to oligotroph and prefers low nutrient availability and low pH 47,48 . Accordingly, we expect that the relative abundance of Acidobacteria increases while Actinobacteria declines in continuous cropping soils. This was not the case in our study. The relative abundance of Sphingobacterii was observed only in crop rotation soils, indicating that it was sensitive to soil environmental disturbances such as N and pH. In terms of both number of OTUs and sequences, the dominant phyla between the abundant and rare sub-community were obviously different (Fig. 2), and more unclassified bacteria were found in the rare taxa. This result was consistent with previous studies in other ecosystems 49,50 . However, Ho et al. 47 found that Copiotrophic Actinobacteria was presented at higher relative abundances in abundant than in rare sub-communities. Nevertheless, this was not the case in our study. This inconsistency may be due to the difference in environmental gradients 31,51 and geography 8,52 . Nowadays, the β-Nearest Taxon Index (βNTI) was broadly used to predict the relative effects of deterministic and stochastic processes on microbial assemblages in microbial ecology 53,54 . Null model analysis revealed that the three sub-communities were driven by stochastic processes in crop rotation soils, while the total and rare microbial taxa were primarily driven by deterministic processes (Fig. 4). Previous results from Jiao et al. 35 and Zheng et al. 55 found that rare microbial taxa were generally affected by deterministic processes in agricultural soils. However, Du et al. 56 thought that rare taxa were more limited by stochastic processes in soil. These inconsistent results may be explained by the different soil environments investigated. We speculate the inconsistency results in our study under different cropping systems may result from the neutral soil pH under crop rotation soils are suitable for most soil microbes, and caused environmental filtering weakened, resulting in the structuring of microbial communities was more influenced by stochastic processes rather than deterministic processes. The variation in environmental factors such as soil N and pH can influence the balance between stochastic and deterministic assembly processes 19 . Wang et al. 57 reported that low N limits the microbial dispersal, and stochastic processes may overwhelm deterministic processes under environmental variation. Future efforts exploring the mechanism of soil environmental factors in driving assemble of soil microbial communities may further proceed.
Compared with the simple description of community structure and diversity indices, network analyses can provide in-depth information on microbial sub-community interactions 14 . In line with our hypothesis, the topological features of abundant and rare taxa in the crop rotation soils were larger than those in continuous cropping soils, except edge, and the results were in line with those of Liu et al. 58 who reported that the co-occurrence network of crop rotation system was more complex than short-term cropping system. This may be due to the plant diversity and the plant residue increase in crop rotation systems, which have an influence on soil microbial communities by altering the exudates secreted by preceding crops 58,59 . This indicates that the crop rotation system has greater co-occurrence networks and more microbial taxa involved in the potential microbial interactions, and can better buffer the environmental changes 60,61 . Jiao et al. 15 reported that abundant rather than rare taxa played vital roles in the construction of microbial networks in oil soils. However, Liu et al. 58 reported that the microbial species that were defined as keystone taxa were not the richest OTUs in the soil under different cropping systems. Nevertheless, in our study, rare taxa in both continuous cropping and crop rotation systems had more connections and closer intra-associations, suggesting rare taxa may exert a great influence on microbial co-occurrences and constitute a more stable community than the abundant biosphere 62 . Overall, these findings indicate that rare species play a unique role in maintaining the function and stability of soil ecosystem.
Environmental factors made different contributions to the assembly of the total, abundant and rare microbial sub-communities in the different cropping soil. The redundancy analysis indicated that within the continuous cropping and crop rotation soils, the total, abundant and rare microbial taxa were well separated (Fig. 3). This finding suggests that cultivation pattern was a main diver in determining the shifts of microbial communities. Previous findings reported by Liu et al. 58 also found that the community of soil microbes was different in response to continuous cropping and crop rotation systems, and the corresponding variations were partially connected with the changes of soil chemical properties. In our study, the variation of bacterial community composition was mainly driven by AP, TOC, SOM, TP, TN, TK, NH 4 + -N, and pH (Fig. 6). Similar conclusions have been drawn by Lauber et al. 63 . and Ai et al. 64 . Our study showed that soil properties made different contributions to the assembly of the rare and abundant microbial taxa in the two different cropping systems ( Table 5). The rare microbial taxa were more connected with soil pH and soil nutrients than the abundant one (Table 5). Soil pH plays the dominant role in determining the soil microbial community structure with different cultivation patterns have been widely reported 51,58 . According to the results reported by Delgado-Baquerizo et al. 65 , the bacterial diversity and composition were partly driven by variation in soil nutrients at a regional scale. This suggests that long-term different agricultural cropping practices can act as a kind of environmental filtering, which can have a strongly impact on soil biotic communities in diverse manners.

Conclusions
Our results showed that the rare and abundant taxa exhibited distinct responses to the different cropping systems. The rare taxa assembly was driven by deterministic processes in continuous cropping soils and by stochastic processes in crop rotation, while the abundant sub-community was driven primarily by stochastic processes in both cropping systems. The α-diversity indices of rare taxa tended to be higher than those of the abundant ones in the two cropping systems. Moreover, the network topologies of rare taxa were more complex than those of abundant taxa in the two cropping systems. Our results suggested that rare taxa functioned as the majority of microbial diversity to enhance bacterial resilience and resistance to continuous cropping disturbances and they served as www.nature.com/scientificreports/ an important role in maintaining the stability in the two different cropping systems. Our study highlights the ecological importance of rare taxa, which can be harnessed in the future for sustainable agriculture production.

Study sites and sampling locations. This research was carried out in Yunnan Academy of Tobacco
Agricultural Sciences' Yanhe Research Farm, Yunnan Province, China (24.14°N, 102.30°E). The study sites are characterized by red soil (classified as Ultisol in USDA (2014) soil taxonomy), which is the dominant type soil for flue-tobacco production in Yunnan 2 . The cultivar for flue-cured tobacco was K326. The mean annual precipitation in this area is about 1160 mm, most between June and October, and the mean annual temperature is 15.9 °C. The altitude in this area is 1680 m, with 2072 h of annual average sunshine. The soil samplings were collected from two tobacco management systems: soil with tobacco monoculture and soil with tobacco-rice rotation. For the continuous cropping soils, the study site was established in 1998, and it has been continuous cropping tobacco for more than 22 years without rotation with other crops. For the noncontinuous cropping soils, the fields were 1 year for tobacco planting and another year for rice cultivation. The two fields were annually added with appropriate chemical NPK fertilizers (97.5 kg ha −1 , N: P 2 O 5 : K 2 O = 1:1:2.5). Other agronomic management measures were same for continuous and non-continuous cropping tobacco fields.
In August 2019, the soil samples were collected from 30 plots (3 × 10 m rectangular plot) after harvesting flue-cured tobacco. Five soil cores from the tillage layer (0-15 cm depth) were collected from each plot using a 5-cm diameter soil corer after removing surface material by hand, and mixed thoroughly and passed through a 2-mm sieve as a single sample. Therefore, there are fifteen soil samples for continuous and non-continuous cropping tobacco fields, respectively. The obtained soil samples were separated into two subsamples: (1) one portion for soil microbial DNA extraction and soil microbial biomass determinations (stored at − 80 °C), (2) the other portion for soil physical and chemical properties determinations. Soil texture, pH, soil organic matter (SOM), total organic carbon (TOC), total nitrogen (TN), total phosphorus (TP), total potassium (TK), dissolved organic nitrogen (DON), nitrate-nitrogen (NO 3 --N), ammonium-nitrogen (NH 4 + -N), available phosphorus (AP) and available potassium (AK), microbial biomass carbon (MBC), microbial biomass nitrogen (MBN), and microbial biomass phosphorus (MBP) were tested using standard methods 66,67 . DNA extraction, 16S rRNA gene sequencing and sequence processing. Briefly, the total soil genomic DNA was extracted from 0.25 g of soil using the PowerSoil® DNA Isolation Kit (MoBio Laboratories, Inc., Carlsbad, USA) according to the manufacturer's protocol. The extracted DNA was stored at − 80 °C until used for Illumina MiSeq sequencing. The bacteria-specific 16S rRNA primers 515F (5′-GTG CCA GCMGCC GCG G-3′) and 907R (5′-CCG TCA ATTCMTTT RAG TTT-3′) were used to amplify the V4-V5 hypervariable regions of bacterial 16S rRNA genes 68 . The amplified PCR products were sequenced on Illumina MiSeq PE300 platform (Illumina Inc., San Diego, CA, USA) at Genesky Biotechnologies Inc., Shanghai, China using 250-bp paired-end reads. The sequencing analysis was as described previously 17,69 . Detailed methods of DNA extraction and purification were available in Jiang et al. 69 .
Raw Illumina paired-end reads were assembled with FLASH 1.2.7 70 , and sequence processing and quality filtering of reads were performed using the pyrosequencing pipeline tools from the QIIME (http:// qiime. sourc eforge. net/) 71 . Reads less than 200 bp and phred score < 25 were discarded, and the high-quality sequences were picked out by the USEARCH using UPARSE software package (http:// drive5. com/ uparse/) 72 . The high-quality sequences with 3% dissimilarity were split into the same operational taxonomic units (OTUs) by the UPARSE pipeline 72 . Subsequently, the representative sequence from each OTU was aligned by PyNAST, and each OTU taxonomic identity was predicted by the ribosomal database project (RDP) classifier using the Greengenes database. Furthermore, the OTUs that contained < 2 reads and singletons were eliminated and all samples were rarefied to 1,244,000 sequences per sample for further analysis. The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2017) in National Genomics Data Center (Nucleic Acids Res 2020), Beijing Institute of Genomics (China National Center for Bioinformation), Chinese Academy of Sciences, under accession number CRA003510 that are publicly accessible at https:// bigd. big. ac. cn/ gsa. Definition of abundant and rare taxa. The abundant and rare OTUs were defined following our previous publications that referred to their local (one sample) and regional (across samples) relative abundances 68 . At the local level (i.e., in one sample), locally abundant operational taxonomic units (OTUs) were defined as the OTUs with relative abundances > 0.1% of the total sequences within a sample, whereas locally rare OTUs were defined as their abundances were < 0.1% within a sample 39 . At the regional level (i.e., across all samples), regionally abundant and rare taxa were defined as OTUs with average relative abundances above and below 0.01%, respectively 17,73 . Finally, the downstream analyses were performed at three levels: whole OTUs, abundant, and rare OTUs.
Sequencing data analysis. Both bacterial α-diversity indices [e.g., the OTU richness and Shannon diversity (H)] and β-diversity based on Bray-Curtis distance between samples were calculated by QIIME to estimate the ecological relationships of bacterial taxa within and among samples, respectively. Furthermore, the same number of sequences from each sample (1,244,000 reads per sample for bacteria) to avoid the effects of different sequencing depths. Nonmetric multidimensional scaling (NMDS) with Bray-Curtis dissimilarities was used to show the differences in the whole, abundant and rare bacterial communities between the continuous and non-continuous cropping tobacco soils using the R vegan package. The microbial community dissimilarities between the continuous and non-continuous cropping tobacco soils were estimated. The distributions of rare www.nature.com/scientificreports/ and abundant bacteria were calculated by Wilcoxon rank-sum tests. Furthermore, the P-values were corrected using the false discovery rate (FDR) correction in all statistical analyses. The microbial phylogenetic assembly on a within-community scale was estimated by the nearest-taxon index (NTI) according to the null model in the Picante package 74,75 . The high or positive NTI values indicate taxa clustering across the overall phylogeny, while low or negative NTI values represent overdispersion of taxa across the phylogeny 76 . On the other hand, the β-NTI index was used to estimate the microbial assembly processes 8 . According to Webb et al. 77 , the β-NTI measures the deviation of the β-mean nearest taxon distance (β-MNTD) from the β-MNTD of the null model, and this was calculated in Phylocom v42. Furthermore, values of |β-NTI| > 2 and |β-NTI| < 2 indicate a community that is dominated by deterministic processes and stochastic processes, respectively 8,53 .
The bacterial co-occurrence patterns of the continuous and non-continuous cropping soils were assessed by network analysis, and networks were visualized according to the interactive Gephi platform (https:// gephi. org/) 22,39 . At the regional level (i.e., across all samples), the co-occurrence networks of rare and abundant bacterial taxa were built based on the correlation analysis for the continuous and rotational cropping soils as described by Jiao et al. 39 . Firstly, we calculated all pair-wise Spearman's correlations between any two OTUs before constructing networks. Secondly, we cut off the two OTUs with Spearman's correlation coefficients less than 0.6 and FDR-corrected p-values more than 0.01 to construct the networks. Each node represents one OUT, and each edge in the network represents significant correlations between two nodes. The node-level topology features were assessed by degree, betweenness, closeness, and eigenvector centrality. The network-level topology features were estimated by a set of indexes: the average path length, network diameter, average degree, graph density, clustering coefficient, and modularity 22 .
The phylogenetic molecular ecological networks of species (pMENs) for the continuous and non-continuous cropping soils were estimated and analyzed by the Molecular Ecological Network Analysis Pipeline (http:// ieg2. ou. edu/ MENA) 78 Detailed information about the network theories, including algorithms, pipeline structure, and operational procedures were described in detail by Deng et al. 78 and Zhou et al. 79 . The major pMENs network properties are as follows: the average geodesic distance (GD), average connectivity (avgK), and average clustering coefficient (avgCC), which were used to test differences among bacterial sub-communities. Besides, the network stability was estimated by the modularity index. The topological roles of nodes were characterized by their positions (within-module connectivity (Zi) and among-module connectivity (Pi)) 68 .
Redundancy discriminatory analysis (RDA) was performed to calculate the influence of the explanatory variables on the bacterial communities following the forward selection procedure 80 . The correlations between bacterial communities (based on Bray-Curtis community dissimilarity) and soil physicochemical characteristics were assessed by Mantel tests.
Approval for plant experiments. Experimental research and field studies on plants (either cultivated or wild), including the collection of plant material, must comply with relevant institutional, national, and international guidelines and legislation.